close all;
clear all;
fclose all;


path1=['C:\Users\Tom�s\Desktop\Temp\'];
% path1=['C:\Users\Tom�s\Desktop\Temp\SR_MLI capacitance traces\25\'];
dir_to_search=[path1]
formatpattern = fullfile(dir_to_search, '/*.ibw');
dInfoNS = dir(formatpattern);
for i=1:numel(dInfoNS)
    nameCell{i}=dInfoNS(i).name;
end
nameCellSort=natsortfiles(nameCell);% natural order sort file names.


%%  make a logical vector for each pattern type so I can perform operations on each group of files, regardless of how they were acquired.
for i=1:numel(dInfoNS)
    temp=IBWread([dInfoNS(i).folder '\' nameCellSort{i}]);
    
    if    numel(strfind(temp.WaveNotes,'Spontaneous'))>0  %if the Wavenotes contain the pattern for each logical
        spontLogical(i)=1;
    else
        spontLogical(i)=0;
    end
    
    if  numel(strfind(temp.WaveNotes,'vStep;a1:-10'))>0   %for capacitance and input resistance
        vStepLogical(i)=1;
    else
        vStepLogical(i)=0;
    end
    
    if   numel(strfind(temp.WaveNotes,'vStep;a1:-30'))>0  % for Ih
        IhLogical(i)=1;
    else
        IhLogical(i)=0;
    end
    
    if numel(strfind(temp.WaveNotes,'IStep'))>0  % for IF curve
        IFLogical(i)=1;
    else
        IFLogical(i)=0;
    end
    
end

IFLogical=logical(IFLogical);
IhLogical=logical(IhLogical);
vStepLogical=logical(vStepLogical);
spontLogical=logical(spontLogical);

%% Analyze capacitance and input resistance
fileNames=nameCellSort(vStepLogical);

for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
    dataMat(:,i)=temp.y;
end

trace= median(dataMat,2);
% this section was taken from the old script, may be useful
Dx=temp.dx; % 1e5 is sampling rate;
DV=-0.01; %10 mV;
t=(0:(numel(trace)-1))*Dx;
figure;
plot(t,trace);
xlabel('Time (s)');
ylabel('I (pA)');
hold on;
plot(0.1,trace(1e4),'r.','markersize',10);
%intStart= find(trace == min(trace(0.8e4:1.2e4))); %the minimum of the trace is the negative peak of the capacitative transient.

[pks,locs]=findpeaks(abs(diff(diff(trace))),'minPeakheight',2); %use the first peak of the second derivative to find onset of cap transient
intStart= locs(1);
hold on;
plot((intStart*Dx),trace(intStart),'r.','markersize',10);
intEnd=1.3e4; %manually selected.
steadyStateIndex=intEnd-500:intEnd+500;
steadyState=mean(trace(steadyStateIndex)); %value to be substracted for integration, only want to integrate the transient and not leak.
baseline= mean(trace(1:1e4));
plot(steadyStateIndex*Dx,trace(steadyStateIndex),'g');
xlim([8e3 intEnd+2000]*Dx);

area(((intStart:intEnd)*Dx),trace(intStart:intEnd),steadyState);

% Create a steady state substracted trace for integration;
SSsub=trace-steadyState;
integral=trapz((1:numel(intStart:intEnd))*Dx,(SSsub(intStart:intEnd)*1e-12));% Trapezoidal numerical integration of the capacitative transient after it has been steady state substracted.
capacitance=round(integral*1e12/DV,1);% capacitance in pF
Rm=round((DV/((steadyState-baseline)*1e-12))/(1e6),1); %Rm in megaohm.
hline(steadyState);
hline(baseline);

title([' Cm=' num2str(capacitance) ' pF, Rm= ' num2str(Rm) ' MOhm' ', Vstep=' num2str(DV*1e3) ' mV']);
obj=scalebar('XLen',0.01,'XUnit','s','YLen',20,'YUnit','pA')

clear fileNames dataMat

%% Analyze spontaneous firing rate (on cell, whole cell)
fileNames=nameCellSort(spontLogical);

for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
   % dataMat(:,i)=medfilt1(temp.y,10);
   
    
figure;plot(temp.y)
prompt = 'set threshold for pk detection (prominence, either negative or positive)         ';
thresh = input(prompt)
close gcf


if i==2
    figure;plot((1:numel(temp.y))*temp.dx*1e3, temp.y,'k')
xlim([0 2000])
obj=scalebar('XLen',100,'XUnit','ms','YLen',20,'YUnit','pA')
axis off
box off
end




% extract spike indices using peakfind
if thresh>0
[pks,locs]=findpeaks(temp.y,'minPeakprominence',thresh,'minpeakDistance',200);
else
[pks,locs]=findpeaks(temp.y*-1,'minPeakprominence',-thresh,'minpeakDistance',200);
end
% calculate mean firing rate and firing rate CV
if numel(locs)>1
FR(i)=1/(temp.dx*(mean(diff(locs))));
CV(i)=std(diff(locs))/mean(diff(locs));

else
    FR(i)=nan;
    CV(i)=nan;
end


% calculate AP mean waveform and half width
spikeIndex=locs(5:end-5);
if numel(spikeIndex)>0
window=5*1e-3/temp.dx;% first term is window in ms
for j=1:numel(spikeIndex)
     mat4STAPre(j,:)= temp.y([(spikeIndex(j)-window):(spikeIndex(j)+window)]);
end

STATracePre(i,:)=nanmean(mat4STAPre,1);
timeSTA=(1:size(STATracePre,2))*temp.dx*1e3;
end
clear mat4STAPre
end
%    
% CV
% FR
if exist('STATracePre')==1
if size(STATracePre,1)==3
%baseline=prctile(STATracePre(3,1:round(window)),5); % use the peak of the
%second derivative instead to find the inflection point.
secDer=diff(diff(STATracePre(3,:)));
[~, SDPkIdx]=max(secDer);
baseline=STATracePre(3,SDPkIdx);
peak=max(STATracePre(3,:));
figure;plot(timeSTA,STATracePre(3,:))
hline(peak,'--k')
hline(baseline,'--k')
halfMax=baseline + 0.5*(peak-baseline);



sub=STATracePre(3,:)-halfMax;
[~,idx1]=min(abs(sub(1:window)))
[~,idx2]=min(abs(sub(window:end)))
idx2=idx2+window;
idx1=idx1*temp.dx*1e3
idx2=idx2*temp.dx*1e3
width=idx2-idx1; % in ms
hold on
 plot([idx1 idx2],[halfMax halfMax],'k')
title({['Width at half max = ' num2str(idx2-idx1) ' ms'],['Firing rate = ' num2str(FR(3)) 'Hz, CV= ' num2str(CV(3))]})
xlabel('time (ms)')
ylabel('Vm (mV)')
else
end
else
    width=nan;
end

clear dataMat

%% Ih section

fileNames=nameCellSort(IhLogical);

for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
    dataMat(:,i)=medfilt1(temp.y,10);
end

figure;plot(t,median(dataMat,2))
ylabel('current (pA)')
xlabel('time (s)')
clear dataMat

%% IF section

fileNames=nameCellSort(IFLogical);

for i=1:numel(fileNames)
    temp=IBWread([dInfoNS(i).folder '\' fileNames{i}]);
    tempNotes=temp.WaveNotes;
    
    idxString=strfind(tempNotes,'a2:');
    
    Step(i)=str2num(tempNotes(idxString+3:end-5)); % need a better way to find the value bc first part can change.
    
    dataMat(:,i)=medfilt1(temp.y,10);
    
    %find spikes and count them.
    [pks,locs]=findpeaks(diff(temp.y),'minPeakprominence',0.1,'minpeakDistance',100);
    
    numSpikes(i)=numel(locs);
    FRSteps(i)=1/(temp.dx*(mean(diff(locs))));
    
end
FRSteps(isnan(FRSteps))=0; %get rid of the nans for plotting

figure;plot(Step,numSpikes,'linewidth',2)
hold on 
plot(Step,numSpikes,'k.','markersize',10)
ylabel('spike number')
xlabel('current injection (pA)')

figure;plot(Step,FRSteps)
hold on 
plot(Step,FRSteps,'ko')
ylabel('firing rate (Hz)')
xlabel('current injection (pA)')





%%simple 1D scatter plot
% 
% figure;
% hold on
% jitt=rand(numel(data),1)/3
% plot(1-1/6+jitt,data,'k.','markersize',30)
% plot([0.5 1.5], [mean(data) mean(data)],'k','linewidth',2)
% xlim([0 2])
% set(gcf,'color','white');
% set(gca,'FontSize', 18);
% box off
% xticks({})

CV
FR
capacitance
Rm


fclose all


